clear;
clc;
close all;
addpath('../02_functions');


%------------------------------------------------------------------------%
% This file creates table 4 in the paper 
%------------------------------------------------------------------------%


% load -----------------------------------------------------------------------

file_name = '../01_data/es_quotes_5minute_EST';

% Call returns routine -----------------------------------------------------------------------

start_date              = datenum(1998,1,5);
finish_date             = datenum(2020,12,29);

winsorizing             = 1;
alpha                   = 0.99;
replace_option          = 1;

results                                                         = make_returns_288(file_name,start_date,finish_date,winsorizing,alpha,replace_option);
log_return_HbH_AD                                               = 100*252*results.log_return_HbH_AD;
log_return_AD                                                   = 100*252*results.log_return_AD;
T                                                               = results.T;
caldt_daily                                                     = results.caldt;

% ----------------------------------------------------------------%
% LOAD TRADES
% ----------------------------------------------------------------%

data = load('../01_Data/es_trades_5minute_EST.mat');

No_trades                                                                   = data.trades_per_min_mat;
No_contracts                                                                = data.trade_quantity_per_min_mat;
Signed_volume                                                               = data.signed_vol_per_min_mat;

beg_TRADES                                                                  = find_closest(data.caldt_daily,start_date);
fin_TRADES                                                                  = find_closest(data.caldt_daily,finish_date);

caldt_trades                                                                = data.caldt_daily(beg_TRADES:fin_TRADES);
No_trades                                                                   = No_trades(beg_TRADES:fin_TRADES,:);
Volume                                                                      = No_contracts(beg_TRADES:fin_TRADES,:);
Signed_volume                                                               = Signed_volume(beg_TRADES:fin_TRADES,:);

% ----------------------------------------------------------------%
% % Load VIX indices
% ----------------------------------------------------------------%

file_name               = '../01_data/VIX_index_5min_EST.mat';
data                    = load(file_name);

beg                     = find_closest(data.caldt_daily,start_date);
fin                     = find_closest(data.caldt_daily,finish_date);

caldt_VIX               = data.caldt_daily(beg:fin);
caldt_VIX_mat           = data.caldt_mat(beg:fin,:);
VIX                     = data.price_mat(beg:fin,:);

indx_beg = find(hour(caldt_VIX_mat(1,:)) == 10 & minute(caldt_VIX_mat(1,:)) == 00);
indx_fin = find(hour(caldt_VIX_mat(1,:)) == 15 & minute(caldt_VIX_mat(1,:)) == 15);
caldt_ID = caldt_VIX_mat(1,indx_beg:indx_fin);
caldt_ID.Format = 'HH:mm';
VIX = VIX(:,indx_beg:indx_fin);

VIX_open = VIX(:,1);
VIX_close = VIX(:,end);

VIX_change1 = VIX_open(2:end) - VIX_close(1:end-1);
VIX_change2 = VIX_close(2:end) - VIX_open(2:end);
VIX_change3 = VIX_close(2:end) - VIX_close(1:end-1);
VIX = VIX_close;
caldt_VIX = caldt_VIX(2:end,1);

% ----------------------------------------------------------------%
% INTERSECT
% ----------------------------------------------------------------%
[C,IA,IB]               = intersect(caldt_trades,caldt_daily);
caldt_daily             = caldt_daily(IB);
log_return_HbH_AD       = log_return_HbH_AD(IB,:);
log_return_AD           = log_return_AD(IB,:);

caldt_trades            = caldt_trades(IA);
No_trades               = No_trades(IA,:);
Volume                  = Volume(IA,:);
Signed_volume           = Signed_volume(IA,:);


[C,IA,IB]               = intersect(caldt_VIX,caldt_daily);
caldt_daily             = caldt_daily(IB);
log_return_HbH_AD       = log_return_HbH_AD(IB,:);
log_return_AD           = log_return_AD(IB,:);
caldt_trades            = caldt_trades(IB);
No_trades               = No_trades(IB,:);
Volume                  = Volume(IB,:);
Signed_volume           = Signed_volume(IB,:);

caldt_VIX               = caldt_VIX(IA,1);
VIX                     = VIX(IA,1);
VIX_change1             = VIX_change1(IA,1);
VIX_change2             = VIX_change2(IA,1);
VIX_change3             = VIX_change3(IA,1);


% ----------------------------------------------------------------%
% Order imbalance calculations
% ----------------------------------------------------------------%

indx                    = [3 : 3 : length(Signed_volume(1,:))];
indx_loop               = [0 , indx];
Signed_volume_15b15_AD  = nan(size(Signed_volume,1),24*4);
Volume_15b15_AD         = nan(size(Volume,1),24*4);

for H = 1 : round(length(Signed_volume(1,:))/12*4)
    Signed_volume_15b15_AD(:,H) = sum(Signed_volume(:,indx_loop(H)+1 : indx_loop(H+1)),2);
    Volume_15b15_AD(:,H)        = sum(Volume(:,indx_loop(H)+1 : indx_loop(H+1)),2);
end

lead = 1;

sum_signedVol_start = 86;
sum_signedVol_end   = 89;

close_signed_volume                      = sum(Signed_volume_15b15_AD(:,sum_signedVol_start:sum_signedVol_end),2);
close_volume                             = sum(Volume_15b15_AD(:,sum_signedVol_start:sum_signedVol_end),2);
indx_no_trading_EOD_volume               = find(close_volume < 100);
RSV                                      = close_signed_volume./close_volume;
Delta_VIX = VIX_change3;


% ----------------------------------------------------------------%
% Loop over positive and negative RSV days seperately
% ----------------------------------------------------------------%

for choose = 1 : 2
    
    if choose == 1
        caldt_tbs       = caldt_daily(1:end-1);
        RSV_tbs         = RSV(1:end-1);
        VIX_tbs         = VIX(1:end-1);
        Delta_VIX_tbs   = Delta_VIX(1:end-1);
        idx_choose      = find(RSV_tbs<0);
        disp('----------------------------------------')
        disp('NEGATIVE RSV DAYS                       ')
        disp('----------------------------------------')
    elseif choose == 2
        caldt_tbs       = caldt_daily(1:end-1);
        RSV_tbs         = RSV(1:end-1);
        VIX_tbs         = VIX(1:end-1);
        Delta_VIX_tbs   = Delta_VIX(1:end-1);
        idx_choose      = find(RSV_tbs>0);
        disp('----------------------------------------')
        disp('POSITIVE RSV DAYS                       ')
        disp('----------------------------------------')
    end
    
    caldt_tbs     = caldt_tbs(idx_choose,1);
    RSV_tbs       = RSV_tbs(idx_choose,1);
    VIX_tbs       = VIX_tbs(idx_choose,1);
    Delta_VIX_tbs = Delta_VIX_tbs(idx_choose,1);
    
    % ----------------------------------------------------------------%
    % Choose sorting variables and set up portfolio sizes
    % ----------------------------------------------------------------%
    
    portfolios = 0 : 33 : 100;
    NN = length(portfolios)-1;
    X1 = [RSV_tbs VIX_tbs];
    
    % ----------------------------------------------------------------%
    % Sorts Number 1
    % ----------------------------------------------------------------%
    
    [caldt_sorted_1 , X1_sorted] = sort_High_Low(caldt_tbs,X1,portfolios);
    
    % ----------------------------------------------------------------%
    % Sorts Number 2
    % ----------------------------------------------------------------%
    
    for i = 1 : NN
        
        temp1_in = caldt_sorted_1{i,1};
        temp2_in = [X1_sorted{i,1}(:,2) , X1_sorted{i,1}(:,1)]; % swap order for second sort
        [temp1_out , temp2_out] = sort_High_Low(temp1_in,temp2_in,portfolios);
        
        for j = 1 : NN
            caldt_sorted{i,j} = temp1_out{j,1};
            X_sorted{i,j}     = [temp2_out{j,1}(:,2) , temp2_out{j,1}(:,1)]; % and swap back
        end
        
        
    end
    
    % ----------------------------------------------------------------%
    % Mean calculations
    % ----------------------------------------------------------------%
    
    for i = 1 : NN
        for j = 1 : NN
            means_X1(i,j) = nanmean(X_sorted{i,j}(:,1));
            means_X2(i,j) = nanmean(X_sorted{i,j}(:,2));
        end
    end
    
    
    asia_start = 1;
    asia_end = 96;
    
    EU_start = 97;
    EU_end = 120;
    
    OD_start = 97;
    OD_end = 108;
    
    CTC_start = 1;
    CTC_end = 267;
    
    TAB_OD = [];
    
    % OD bucket mean calculations
    
    for i = 1 : NN
        for j = 1 : NN
            caldt_ij      = caldt_sorted{i,j};
            [C,IA,IB]     = intersect(caldt_daily,caldt_ij);
            idx           = IA+1;
            TAB_OD(i,j)   = nanmean(nansum(log_return_AD(idx,OD_start:OD_end),2));
        end
    end
    
    % High minus Low and p-values
    
    row_stats = [];
    
    % differences across rows
    for i = 1 : NN
        
            caldt_i1      = caldt_sorted{i,1};
            [C,IA,IB]     = intersect(caldt_daily,caldt_i1);
            idx           = IA+1;
            temp1         = nansum(log_return_AD(idx,OD_start:OD_end),2);
            
            caldt_i3      = caldt_sorted{i,3};
            [C,IA,IB]     = intersect(caldt_daily,caldt_i3);
            idx           = IA+1;
            temp2         = nansum(log_return_AD(idx,OD_start:OD_end),2);
            
            [~,p] = ttest2(temp2,temp1);
            row_stats(i,1) = nanmean(temp2)-nanmean(temp1);
            row_stats(i,2) = p;
    
    end
    
    col_stats = [];
    
    % differences across columns
    for i = 1 : NN
            caldt_1i      = caldt_sorted{1,i};
            [C,IA,IB]     = intersect(caldt_daily,caldt_1i);
            idx           = IA+1;
            temp1         = nansum(log_return_AD(idx,OD_start:OD_end),2);
            
            caldt_3i      = caldt_sorted{3,i};
            [C,IA,IB]     = intersect(caldt_daily,caldt_3i);
            idx           = IA+1;
            temp2         = nansum(log_return_AD(idx,OD_start:OD_end),2);
            
            [~,p] = ttest2(temp1,temp2);
            col_stats(1,i) = nanmean(temp1) - nanmean(temp2);
            col_stats(2,i) = p;
    
    end
    
    
    
    
    if choose == 1
        TAB_OD   = flipud(TAB_OD);   % to have the numbers in the intuitive direction
        means_X1 = flipud(means_X1); 
        row_stats = flipud(row_stats); 
    end
    
    means = [100*means_X1 , means_X2];
    OD = NaN(5,6);
    OD(1:3,1:3) = TAB_OD;
    OD(4:5,1:3) = col_stats;
    OD(1:3,4:5) = row_stats;
    
    % ----------------------------------------------------------------%
    % Print Results
    % ----------------------------------------------------------------%
    
    x_string = {'VIX Low' , 'VIX Med' 'VIX High' , 'VIX Low' , 'VIX Med' 'VIX High'};
    y_string = {' ' , 'RSV Low' , 'RSV Med' 'RSV High'};
    

    disp('RSV bucket means ... VIX bucket means:')
    info.cnames = strvcat(x_string);
    info.rnames = strvcat(y_string);
    info.fmt = '%6.2f';
    info.swidth = 220;
    info.hspc = 5;
    mprint1(means, info);

    
    x_string = {'VIX Low' , 'VIX Med' 'VIX High','High-Low','p-value',' ' };
    y_string = {' ' , 'RSV Low' , 'RSV Med' 'RSV High','High-Low','p-value'};

    disp('OD bucket means:')

    info.cnames = strvcat(x_string);
    info.rnames = strvcat(y_string);
    info.fmt = '%6.2f';
    info.swidth = 220;
    info.hspc = 5;
    mprint1(OD, info);

    
    
    
end